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The stability of sand castles is determined by the structure of wet granulates. Experimental data 
about the size distribution of fluid pockets are ambiguous about their origin. We discovered that 
contact angle hysteresis plays a fundamental role in the equilibrium distribution of bridge volumes, 
and not geometrical disorder as commonly conjectured, which has substantial consequences on 
the mechanical properties of wet granular beds, including a history dependent rheology and lowered 
strength. Our findings are obtained using a novel model where the Laplace pressures, bridge volumes 
and contact angles are dynamical variables associated to the contact points. While accounting for 
contact line pinning, we track the temporal evolution of each bridge. We observe a cross-over to a 
power-law decay of the variance of capillary pressures at late times and a saturation of the variance 
of bridge volumes to a finite value connected to contact line pinning. Large scale simulations of 
liquid transport in the bridge network reveal that the equilibration dynamics at early times is well 
described by a mean field model. The spread of final bridge volumes can be directly related to the 
magnitude of contact angle hysteresis. 
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I. INTRODUCTION 

Our daily life is strongly affected by the mechanical 
properties of wet granulates, be it through the stability 
of the soil on which we construct our buildings or the 
consistency of food we enjoy. The structure of the fluid 
interfaces in a wet granular bed essentially determines the 
strength of capillary cohesion and thus plays a key role for 
the prediction and mitigation of devastating events such 
as landslides or the clogging of industrial pipes 00 - 
For a thorough understanding, a realistic model of fluid 
movement is indispensable: The capillary pressure driven 
fluid transport in wet granulates 0 H changes the fluid 
distribution which in turn suggests that the mechanical 
properties are affected by the fluid displacement. 

Indeed, shear thinning for example is observed for par¬ 
tially saturated beds of glass beads 0 , 0 , 0 due to a 
delayed reconfiguration of fluid interfaces. Knowledge 
about the evolution of the fluid distribution in the gran¬ 
ular assembly is relevant to many industrial processes 
that involve coating of powders and grains Q, or mixing 
additives with particulate materials 01]. The precise 
microscopic features which crucially determine the na¬ 
ture of the equilibration and fluid transport processes 
are yet largely unexplored. 

In this article we discover in simulations of the fluid 
equilibration process in a bed of spherical beads that 
contact line pinning can have a profound effect on the 
evolution of capillary bridge volumes. In agreement with 
the majority experiments on fluid transport and pressure 
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equilibration in wet granular beds a i i na , we con¬ 
sider the wetting fluid to be the minority phase forming 
capillary bridges at grain contacts or in narrow gaps be¬ 
tween neighboring grains. Our simulations demonstrate 
that contact angle hysteresis broadens the distribution 
of final bridge volumes while the polydispersity of bead 
radii and the distribution of gap separations in the granu¬ 
lar bed appear to play only a minor role. The present ap¬ 
proach accounts for the discreteness of the bridge network 
and incorporates the history dependence of the bridge 
shape caused by contact angle hysteresis. 

Combining the model for fluid transport with the Con¬ 
tact Dynamics (CD) algorithm [lllTld] ] allows us to sim¬ 
ulate the motion of beads due to an external driving and 
the redistribution of liquid after the external energy in¬ 
put has been stopped. In this way, a direct comparison of 
our numerical results to the evolution of bridge volumes 
using X-ray microtomography imaging of fluid distribu¬ 
tion measured by Scheel et al. 0 , 110 ] is possible. In these 
experiments vertical shaking of the granular bed was em¬ 
ployed to distribute the wetting liquid on the surface of 
the beads. Shortly after the agitation of the wet granu¬ 
lar bed (e.g. by stirring of shaking) has been stopped, the 
surfaces of the beads are covered with a ‘thick’ macro¬ 
scopic film of the wetting fluid 0 ]. Frequent closing and 
opening of contacts leads to a uniform distribution on the 
beads since the timescale of collisions is typically much 
smaller than timescales of viscous flows on the surface 
of the beads or of diffusion through the continuous fluid 
phase. 

After the beads have settled into stable positions in the 
bead pack, the wetting fluid starts of flow to the points 
of closest proximity of opposing bead surfaces, forming 
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FIG. 1 . (Color online) Effect of contact angle hysteresis on 
the equilibration dynamics of two capillary bridges labeled i = 
1, 2 between beads with identical radii in mechanical contact. 


macroscopic capillary bridges. At the end of this process, 
the liquid volume varies from bridge to bridge reflecting 
both the disorder in the bead pack and the distribution 
of fluid patches on the beads immediately after shaking. 

The initial regime of bridge formation is followed by an 
equilibration of capillary pressure which must involve an 
exchange of wetting fluid between neighboring bridges. 
In this ‘equilibration regime’ the majority of the wetting 
fluid will be contained in macroscopic capillary bridges, 
cf. the sketch in Fig. [T] In the present work we will fo¬ 
cus on the collective evolution of bridge volumes in the 
equilibration regime where we assume that volume con¬ 
servation holds for the wetting fluid. Our model for the 
equilibration dynamics is based on the assumption that 
the fluid transport is entirely driven by differences in the 
local capillary pressure. For wetting liquids with a low 
vapor pressure this transport proceeds mainly through 
thin wetting films in the roughness of the beads d, @j. 
A diffusive transport through the vapor phase has to be 
considered for volatile liquids that exhibit a sufficiently 
high vapor pressure [ijl 0 • 

Central aim of the present work is to quantify the im¬ 
pact of contact angle hysteresis and polydispersity of the 
particles in the granular bed onto the equilibration of 
bridge volumes. In particular, we address the question 
for typical timescales of the transport processes. These 
timescales of capillary pressure equilibration play an im¬ 
portant role in the rheology of partially saturated wet 
granular beds. 


Throughout the following considerations, we will ne¬ 
glect contributions of the hydrostatic pressure to the 
pressure difference P = pi — p v between the liquid ( 1 ) 
and vapor (v). Hence, the capillary pressure P is the 
same in every point of the liquid-vapor interface. Pro¬ 
vided a uniform interfacial tension 7, the Young-Laplace 
equation 


P = 2H7 


(1) 


implies that the interface is a surface of constant mean 
curvature H. The mean curvature H is the sum of the 
two principal curvatures n\ and «2 in a respective point 
of the interface. 

The influence of gravity or buoyancy on the capillary 
pressure P will become apparent whenever the maximal 
difference in hydrostatic pressure is of the order or ex¬ 
ceeds 7/-R0, where Rq is the average bead radius. This is 
the case if the vertical extension of the system becomes 
comparable, or exceeds the length 


L± 


7 
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( 2 ) 


where A p = pi — p v is the difference in mass density 
between the wetting liquid l and the ambient fluid phase 
v, while g is the acceleration of gravity. 


A. Liquid transport 

Two different models for the local transport of wetting 
liquid between neighboring capillary bridges can be con¬ 
sidered. In case the wetting liquids forms a thin film on 
the surface of the grains, any pair of bridges on the same 
bead are hydraulic connected. Assuming a quasi-steady 
flow in these films the total flux between two capillary 
bridges on the same grain is proportional to the differ¬ 
ence in their capillary pressure. The corresponding liquid 
mobility, or conductance coefficient depends on the ge¬ 
ometry and distance of the three phase contact lines of 
the bridges. Besides the advective transport mechanism, 
we have to account for a capillary pressure equilibration 
by mass diffusion of the wetting liquid through the am¬ 
bient continuous phase. The latter transport mechanism 
may be effective for wetting liquids with a high vapor 
pressure or liquids with a partial miscibility in the con¬ 
tinuous liquid phase. 


II. PHYSICAL MODEL 

For simplicity, we consider the wetting liquid that 
forms the capillary bridges (minority phase) to be sur¬ 
rounded by a continuous vapor phase (majority phase). 
Analogous statements can be obtained for binary mix¬ 
tures of partially miscible liquids where the minority 
phase is a wetting liquid while the majority phase is a 
non-wetting liquid. 


1 . Viscous film flow 

For systems where the non-volatile liquid forms a wet¬ 
ting film on the surface of the beads we will regard the 
average thickness hg of the film to be insensitive to the 
pressure difference P between the wetting liquid and the 
vapor phase. This assumption is justified whenever the 
typical horizontal length scales of the roughness is much 
smaller than the radius of curvature of the menisci in the 
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macroscopic capillary bridges. In addition, we need to 
assume a sufficiently small microscopic contact angle of 
the liquid-fluid interface on the rough grain surfaces that 
permits the formation of a percolating liquid film [M2. 

Assuming a stationary viscous flow in the thin film at 
any instance in time, the volume flux Qjj between a pair 
of bridges i and j on the same beads will be proportional 
to the difference of capillary pressures Pi — Pj and the 
liquid mobility ~ hp/ij in the thin film fl8l ]: 

Qij = (Pi - Pj) , (3) 

where 77 is the dynamic viscosity of the wetting liquid, 
ho the thickness of the wetting film, and Pij the Laplace 
pressure of the two neighboring bridges. Here, we assume 
that pressure gradients occur over a typical distance of 
the order ~ Rq. Similarly, we estimate the circumference 
of the three phase contact line of the bridges on a bead to 
be of the order ~ Rq. The dependence on the particular 
geometry of the neighboring capillary bridges is adsorbed 
in the dimensionless conductance coefficient Cij , being of 
the order of unity. 

An exact computation of the conductance coefficients 
requires solutions to the Laplace equation of the local 
capillary pressure in the film. Dirichlet boundary data 
for the capillary pressure are given on the contact lines. 
Conformal mapping of the surface of the sphere to a plane 
allows us to reduce the problem to the computation of 
the capacitance per length of two parallel and infinite 
cylindrical conductors [ 22 j. This results shows that the 
dependence on the size and distance of the bridges is at 
most logarithmic. 


2. Diffusion through ambient fluid 


A rough estimate of the diffusive flux of a volatile wet¬ 
ting liquid through its ambient vapor phase can be de¬ 
rived from Fick’s first law E 3 E]. In full analogy to the 
second regime of capillary equilibration in the thin film 
model described above, we will assume that the concen¬ 
tration profile of liquid molecules in the vapor phase has 
already leveled in every single pore of the bead pack. 
Gradients in the concentration are noticeable only on 
length scales much larger than the typical dimension of a 
pore. As shown in the appendix IVl we can write Kelvin’s 
equation in the form 


P n. 


*2 


n v - n„ = 


p*nj 


(4) 


which allows us to express the particle density of liquid 
molecules n v in the vapor phase through the vapor pres¬ 
sure p* and the densities n* and n* at bulk coexistence, 
and the pressure difference P = pi — p v across the curved 
interface. As shown in the appendix El Kelvin’s equa¬ 
tion 0 still holds if we replace the vapor pressure p* and 
densities n*, n* at bulk coexistence of the pure phases 


by the partial vapor pressure p v , and the particle densi¬ 
ties n v , hi for a given atmospheric pressure po of the gas 
phase, respectively. 

Using Fick’s first law, we can express the volume flux 
of liquid caused by diffusion of liquid molecules between 
bridge i and j as 


_ ^ DnfRp 
u u p*n * 2 , 


(Pi-Pj), 


(5) 


where D is the diffusion constant of molecules in the 
vapor phase, cf. El The dimensionless prefactor Cij ac¬ 
counts for the specific geometry of the pore space and 
the bridge interfaces. 

In full analogy to the dimensionless conductance coeffi¬ 
cient Cij in the viscous film model, we expect the prefac¬ 
tor Cij to be of order unity. A computation of the matrix 
elements Cij involves solutions of the three dimensional 
Laplace equations for the stationary density profiles with 
Dirichlet boundary data on the liquid-vapor interface of 
the capillary bridges. 


3. Transport timescales 


In view of the numerical simulations, it is useful to 
non-dimensionalize all relevant physical quantities. In 
the following we employ the average bead radius Rq as a 
unit of length. After rescaling the pressure by 'y/Rp, the 
relation 0 for the volume flux yields a capillary-viscous 
time scale of the form 


= 

~ iK 


( 6 ) 


for the pressure equilibration by fluid transport through 
a thin viscous film. Employing the relation 0 instead 
of 0 leads to a time scale 


_ P*nfRl 

infD 


(7) 


for a diffusive transport through the vapor phase. 

Volume fluxes of both transport modes can be super¬ 
imposed as long as the linear relations between the re¬ 
spective volume flux and pressure difference eqns. 0 and 
© are applicable. This implies that the timescale 


T v T d 
T v + T d 


( 8 ) 


accounts both for flows through thin liquid films and for 
diffusion through the vapor phase. The corresponding 
dimensionless characteristic number 


TV _ D 77 n* 2 Rp 
Tfi p* n* 2 hg 


allows us to discriminate between two different transport 
regimes: for J 1 , the diffusive transport through the 
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vapor phase dominates while for J < 1 the main trans¬ 
port is by viscous flows through the thin wetting film. 
It has been shown experimentally that both cases may 
indeed be encountered fl9l - [2H . The central message here 
is that although the relative importance of the two mech¬ 
anisms may vary and can be expressed by the value of J, 
the form of the transport equations remains unchanged. 
Hence it is sufficient to treat the film flow case, as we will 
do in what follows. 


where the dimensionless factor £ = R c /Rq accounts for 
the Derjaguin mean 


_ 2R1R2 
c = R 1 +R 2 


of bead radii [24| . 


(14) 


C. Contact angle hysteresis 


B. Capillary bridges 

As in a previous Letter [l4{ the capillary pressure P 
of a bridge with volume V spanning a gap with separa¬ 
tion S between beads of equal radius Rq is interpolated 
from tabulated values. To account for contact angle hys¬ 
teresis, we consider capillary bridges with either pinned 
or freely moving contact lines. A suitable parameter to 
describe the position of the pinned contact line of a cap¬ 
illary bridge is the opening angle /3, cf. also the sketch in 
Fig. [0 

The capillary pressure of a bridge is a function 
Po(V,S,P ) for the bridges with pinned contact lines de¬ 
pending, besides on the volume V and gap separation S, 
on the opening angle f3. For bridges with a freely sliding 
contact line and a prescribed contact angle 9, the cap¬ 
illary pressure is expressed as a function P 0 (V,S, 8 ). In 
addition to the capillary pressure, we compute the con¬ 
tact angle 9(V,S, (3) of bridges with a pinned contact line 
and the opening angle (3 (V,S, 8 ) for those with a slid¬ 
ing contact line. To numerically compute the functions, 
Po(V,S,0), P 0 (V,S, 8 ), 6 {V,S,P), and /3(V,S,9) we em¬ 
ployed numerical energy minimizations using the public 
domain software Surface Evolver [23| where we consid¬ 
ered only axially symmetric interfacial profiles using an 
effectively two dimensional representation of the bridge 
state. The rupture distance of a bridge with fixed volume 
V is computed from the approximate expression 

-So* = (i + 0 V 1/3 (10) 

as given by Willct et al. [24| . As shown in one of the fol¬ 
lowing works, Ref. [aJ. eqn. JED is still a good approx¬ 
imation for bridges with a pinned contact line provided 
that the opening angle f3 is not too small. In this case, 
the actual contact angle 9 = 0(V, S, j3) of the bridge is 
used in eqn. JED- 

To account for the dependence of capillary pressure P 
on the volume V, the surface to surface separation S , and 
on the individual radii Ri of the two beads i=l, 2 we use 
scaled quantities 

p(v,s,9) = r 1 p 0 (r 3 v,r 1 s,9) (n) 

p{y,s,p) = r 1 Mr 3 v,r 1 s,p) m 

S*(V 1 9)=ZS*(C 3 V,9) (13) 


Roughness or variations in the chemical composition 
of the surface lead to a history dependent contact angle 
fl8j ]. Assuming ideal surfaces with uniform surface het¬ 
erogeneities the contact angle 9 of a fluid interface in a 
mechanical equilibrium falls into an interval [ 9 r , 8 a ] lim¬ 
ited by the receding and advancing contact angle 9 r and 
9 a , respectively. The contact line starts to recede once 
the local contact angle 9 equals 9 r . Likewise, the contact 
line starts to advance if 9 equals 9 a For all intermedi¬ 
ate values of 9 € [ 9 r . 9 a ] the contact line is immobilized, 
i.e. ‘pinned’. The magnitude of contact angle hysteresis 
is defined as the width A 9 = 9 a — 9 r of the range of static 
contact angles. 

The sketch in Fig. |T] illustrates the impact of contact 
angle hysteresis on the equilibration dynamics for the 
example of two capillary bridges i = 1,2. Starting from 
an initial state with volumes V® > V^°, opening angles 
f3i > (3% and contact angles 0° = 9 a , 9% = 9 r at t = 0, 
the capillary pressure difference P[ J — P 2 >0 drives a 
fluid flux Q 12 > 0 through the thin film between bridge 
1 and 2. As the total volume of bridge 1 and 2 must be 
conserved we have V\ = — V 2 = —Q 12 , using the dot as 
a short hand notation for the total time derivative. Al¬ 
though the bridge 1 is shrinking and bridge 2 is growing, 
all contact lines remain first in a pinned state {f3\ = /3°, 
/?2 = (3 2 i and 81,2 6 [ 9 r , 0 a ]). Only while approaching the 
final state at t —> 00 with asymptotically equal pressures 
P£° = P£°, the contact lines become depinned (9^° = 0 r , 
9^° = 9 a ), but still with bridge volumes Vj_°° > V£°. 
In the absence of contact hysteresis, one would observe 
equal final volumes Vf° = V 2 °°. 

D. Network model 

The dynamics of fluid transport in the static bed of 
beads is studied in a coarse grained network model where 
each capillary bridge forms a node i of the network. The 
nodes are connected by bonds representing the thin 
film between bridges j and i on the same bead. Mass 
conservation of the wetting fluid demands that the rate 
of volume change of bridge i is the sum 

Vi = - J2 Qv= E CijiPj-Pi), ( 15 ) 

over the set M(i) of all bridges j that are connected to 
bridge i. The relative geometry of a pair of neighboring 
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bridges (i,j), the effect of other bridges, or individual 
bead radii is collected in the dimensionless conductance 
coefficient Cij, which is of the order of one. The loga¬ 
rithmic dependence of Cij on the minimum distance of 
contact lines applies only to capillary bridges close to co¬ 
alescence. Since we exclude this exceptional case, we set 
Cij = 1 throughout. 

For given derivatives Vi of bridge volumes eqn. m 
and gap separations Si, we are able to compute the total 
time derivatives of the opening angle Pi and contact angle 
Oi of bridge i. Using the differentials 

a = d v /3(Vi,Si,6i) V + d s P(Vi, Si, Oi) A 
b = d v e(Vi, Si, pi) Vi + d s 9(Vi , a, A) A 


we can express the time derivatives Oi and Pi as 


A = a , 

Pi = 0, 


Sh 

<2 

O 

II 

p p 

IA IV 
o o 

and 0i = 0 a 
and 0i = 0 r 

(16) 

f—i 
£ 

iO 

II 

•qT 

0 r <c 0i <c 0 a . 

(17) 


The dependence on the gap separation S in the differen¬ 
tials becomes important if the beads are allowed to move 
relative to each other. In a static packing, however, the 
angles A and A depend on V, only. 

According to eq. (fTBl) . the contact line of a bridge i 
can only slide outward (A > 0 ) for a contact angle Oi = 
0 a . Capillary bridges with inward sliding contact lines 
(A < 0 ), however, have to satisfy Oi = 0 r . Equations 
m guarantees that the contact angle Oi cannot become 
larger than 0 a or smaller than 0 r . The complementary 
case of an immobilized contact line is taken into account 
by eqs. (ED: Any intermediate value 0 r < 0 t < 0 a of the 
contact angle Oi implies a pinned contact line (A = 0 ). 

As a consequence of the dynamics described by 
eqs. m and ED, the state of a capillary bridge is not 
uniquely determined by the volume V and separation 
S. Depending on whether the contact line is advancing, 
pinned, or receding we have to employ different expres¬ 
sions for the capillary pressure: 

r P{Vi,Si,0 a ) for A > 0 
Pi = \ P(Vi,Si,pi) for P = 0 . (18) 

[ P{V i ,S l ,0 r ) for A <0 

The evolution of liquid volume Vpt), opening angle ppt), 
and contact angle Opt) for each individual bridge i is 
obtained from an integration of eqns. (EDDIED- 

The dependence of the capillary pressure P on the his¬ 
tory of bridge volume V for a fixed separation S is illus¬ 
trated in Fig. [2] The two flat branches in the sketch in 
Fig. [5] correspond to the case of an advancing (solid line 
top) and receding contact line (solid line bottom), while 
the step curves (dashed lines) describe the pressure of 
a capillary bridge with pinned contact lines. Owing to 
the history dependent dynamics specified in eq. (UHP , a 
capillary bridge moving on the branch corresponding to 
an advancing contact line can be encountered only for 



FIG. 2. (Color online) Schematic of the capillary pressure 
P against the volume V of a bridge. The lower and upper 
branch (dashed lines) corresponds to a receding (lower) and 
advancing (upper) contact line. The left and right branches 
(solid lines) describe bridges with a pinned contact line with 
small (left) and at large (right) opening angle. 


growing volumes. Similarly, bridges on the lower branch 
corresponding to a receding contact line are encountered 
only for shrinking volumes. In contrast, the volume of a 
bridge with a pinned contact line can be either growing 
or shrinking. 

A substantial simplification to the full network model 
eqs. (fl5l) can be obtained by replacing the capillary pres¬ 
sure Pj in eqn. m with j € A f(i) by the mean ( P) 
taken over all bridges in the network. The temporal evo¬ 
lution of individual bridge volumes V) in this ‘mean field’ 
approximation is 


Vi = C (N C )((P) - Pi) (19) 

where C is a global conductance coefficient and (N c ) the 
average bridge coordination on a bead. 


E. Simulations 

In order to reproduce the experimental conditions of 
the equilibration experiments in Refs. ins , we simu¬ 
late the particle motion during the preparation step using 
the contact dynamics (CD) algorithm. Since the beads 
employed in the experiments are a sieving fraction of a 
wide size distribution it is reasonable to assume bead 
radii which are uniformly distributed between a mini¬ 
mum radius, R-, and a maximum radius, i? + . The 
degree of polydispersity is parameterized by the ratio 
r po i = R-/R+ < 1, while is employed as our unit 
of length, Rq. The open simulation box is bounded by 
a rigid movable wall at the bottom z = 0 , applying pe¬ 
riodic boundary conditions along the x and y directions. 
Gravity acts on the particles into — z direction with an 
acceleration of gravity g = 9.81m/s 2 . 
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FIG. 3. (Color online) Relaxation of individual capillary 
bridge volumes as a function of time in experiments (symbols) 
taken from Ref. M and in our simulations for Ad = 25° in 
the full network model (‘NW’, black dashed lines) and mean 
field model (‘MF’, red solid lines). The inset shows the equi¬ 
libration dynamics in the network model with Ad = 0°. The 
receding contact angle is d v = 7° in both cases. 


To create a static packing similar to those in the time 
resolved x-ray tomography experiments in Refs. B0, 
the bed of beads is fluidized by an oscillatory motion of 
the lower wall of the container. The wetting fluid is re¬ 
distributed exclusively via rupture and formation of cap¬ 
illary bridges. If two beads touch each other, a small 
bridge is created with 9 = 6 a , where the required, mini¬ 
mum fluid volume is provided by all neighboring bridges. 
Similarly, if the bead separation exceeds the critical dis¬ 
tance S*(V), the bridge ruptures and the volume V is 
equally split and distributed onto all bridges on the two 
beads. Capillary forces are small compared to the typical 
bead mass and acceleration amplitudes during fluidiza¬ 
tion and are therefore neglected. 


III. RESULTS AND DISCUSSION 

The simulated evolution of bridge volumes in a static 
bed of beads after shaking (set up A) is shown in Fig. [3] 
Following the analysis of the x-ray tomography data in 
Ref. 0 , we classify each bridge according to their binned 
volume at time t = 0. Volume averages over all bridges 
being in the same bin at t = 0 are plotted against the time 
t elapsed after the agitation has been stopped. To match 
the parameter of the experiments reported in Ref. 0 , 
we set the polydispersity to r po i = 0.8 and chose a total 
fluid volume corresponding to an overall fluid content of 
W = 2 x 10 -2 of the wetting fluid with respect to the 
total sample volume. The receding contact angle in our 
simulations is set to 9 r = 7° throughout, being the small¬ 
est value that allows for a reliable interpolation between 
the tabled functions P(S,/3,9) and V(S,f3,9). The dis¬ 
tribution of final bridges volumes for an advancing con¬ 
tact angle 9 a = 32° (corresponding to a contact angle 
hysteresis of A 9 = 25°) provides the best match to the 
experimental data in Ref. 0- 

The plot in Fig. [3] demonstrate that the average bridge 
volume of a bin does not need to be a monotonously in¬ 
creasing or decreasing function of time. In particular, 
we observe that some capillary bridges with intermedi¬ 
ate volume first shrink for a short time and then tend 
to grow again. The inset of Fig. [3] shows the case of 
a vanishing contact angle hysteresis A 9 = 0 °, i.e. for 
contact angles 9 a = 9 r = 7°. For asymptotically large 
times all bridge volumes V) virtually converge to the same 
value ( V) = V°°. This indicates that local the geometry 
around contact points and small gaps, i.e. the possible lo¬ 
cations of capillary bridges, are very similar. Because the 
capillary pressure will finally be identical in all bridges, 
the final value V°° is determined only by the contact an¬ 
gle, the number of contacts, and the total liquid volume 
in the granular bed. 


A. Transport timescales 


Once the agitation is stopped we let the beads settle 
by gravity into a mechanically stable packing (setup A). 
The distribution of fluid volumes obtained with this pro¬ 
cedure is close to a decreasing exponential function. To 
extend our calculations to larger systems while saving 
the computational costs for the initial preparation, we 
considered a three dimensional cubic lattice of beads in 
contact applying periodic boundary conditions to the box 
with linear dimensions L = 200 (setup B). To mimic the 
initial state of bridges after the preparation step in A, 
we initialized the bridges with an exponentially decay¬ 
ing probability distribution function (PDF) of volumes 
oc exp(— V/{V)) with a given average volume (V) and 
set all contact angles 9 = 9 a - To follow the evolution 
of individual bridge states, we integrated eqs. (I15II18I) for 
setup A and B numerically in time using a simple forward 
Euler time scheme. 


In addition to the values the contact angles 9 a and 9 r , 
we can employ the time scale To as a free parameter to 
fit the results of our numerical simulations to the exper¬ 
imental data. The best fit between the data of the equi¬ 
libration experiments in Ref. 0 and the simulations in 
set up A shown in Fig. [3] is obtained for To = 2.1 x 10 5 s. 
Rather unexpected, the volume equilibration in set up A 
appears to be almost completed already after t > 10~ 3 
expressed in the dimensionless units of our simulations. 
Rescaling this time by To leads to ~ 250 s for the cross¬ 
over which is corroborated by visual inspection of Fig. [3] 
If we assume that the transport proceeds exclusively 
through viscous flows in thin films, we can employ 
eqn. © to estimate the equivalent thickness of the thin 
film. Setting the bead radius to Rq = 280 /im, the in¬ 
terfacial tension and dynamic viscosity of the aqueous 
Znl 2 solution to 7 = 72mNm _1 and r] = 10 -3 Pa s, re- 
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B. Final volume distribution 
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FIG. 4. (Color online) Probability distribution function 
(PDF) of final bridge volumes V for vanishing contact an¬ 
gle hysteresis according to set up A. Contact angle 9 = 7°. 
Inset: Variance of volumes in the final state for different poly- 
dispersities r po i. 


spectively, we obtain an estimate of ho = 74 nm for the 
effective thickness of the thin wetting film. This values 
agrees well with the measured RMS roughness of glass 
beads by Utermann et al. [26] and the value of ho pro¬ 
posed by Lukyanov et al. for their system of the non¬ 
volatile liquid Tris(2-ethylhexyl) phosphate in a pack or 
quartz grains (Ottawa sand) 0. 

An estimate of the equilibration timescale To in pres¬ 
ence of a purely diffusive transport in the vapor phase 
can be obtained from eqn. ©. Assuming pure water 
as a wetting liquid and ambient conditions with a tem¬ 
perature of 300 K, we set the vapor pressure to p* = 
3.52 x 10 3 Pa, the molar densities to n* = 5.56 x 10 4 mol 
m -3 and n* = 1.41 mol m -3 , the diffusion constant 
to D = 2.57 x 10 _ 5 m 2 s _1 , and obtain a timescale of 
Td = 6.5 x 10 7 s. The low value of the characteristic 
number J = 3.23 x 10 -3 according to eqn. © suggests 
that a transport by thin film flow is favored over a diffu¬ 
sive transport in the vapor phase. A further quantitative 
discussion requires the numerical computation of the ma¬ 
trix elements Cij and Cij , and is left to future work. 

A further argument that supports viscous flows as the 
dominant transport mode in Ref. 0 is due to the high 
molarity of the aqueous Z 11 I 2 solution. Any diffusive 
exchange of only the volatile solvent between bridges, 
excluding and exchange of the non-volatile ionic solute 
leads to gradients in chemical potential of the solvent 
molecules. A simple estimate using Raoult’s law demon¬ 
strates that any gradient in the chemical potential due 
to differences in capillary pressure is small compared to 
the gradients caused by differences in molarities of the 
solute. Hence, a diffusive exchange of bridge volumes is 
largely suppressed for high salt concentrations. 


An explanation of the width of final bridge volumes 
by hydrostatic pressure can be excluded for the equili¬ 
bration experiments by Scheel et al. 0 - Given the in¬ 
terfacial tension 7 and density difference Ap of the 1 M 
aqueous Znl 2 solution to the surrounding air, we compute 
from eqn. © a maximal vertical extension of L j_ ss 1 cm. 
Since the analyzed field of view in the x-ray tomography 
imaging had a height of less than t « 5mm 0 , we con¬ 
clude that only differences in the Laplace pressure were 
driving the exchange of wetting fluid between neighbor¬ 
ing bridges. Gradients of the final distribution of bridge 
volumes caused by a mismatch in the hydrostatic pres¬ 
sure in the bulk phases will contribute to the width of 
the volume distribution only if t > L± holds. 

A remaining finite width of the final volume distribu¬ 
tion can be related to the degree of polydispersity r po i 
of bead radii. Figure [I] shows the probability density 
functions of bridge volumes V obtained in our simula¬ 
tions for a number of size polydispersities r po i at late 
times. Apparently, all PDFs approach a characteristic 
triangular shape while their width is increasing as r po i 
becomes smaller. As expected, the variance of bridge 
volumes, a v = (( V 2 ) — (P ) 2 ) 1 / 2 that quantify the width 
of the PDFs, increases monotonously with r po i, the de¬ 
gree of polydispersity (cf. the inset of Fig. [3j. This weak 
sensitivity is a clear indication that the final distribu¬ 
tions of volumes observed in the experiments cannot be 
explained by the polydispersity of bead radii alone. Cap¬ 
illary bridges between spheres not in mechanical contact 
are rare in packings prepared by setup A and the PDF 
of gap separations S between neighboring beads has a 
negligible impact on the final volume distribution. 

An obvious source for the finite width of the final distri¬ 
bution of bridge volumes is a history dependent contact 
angle. Figure [5ji) shows the PDF of bridge volumes at 
different times after the agitation has been stopped for a 
finite contact angle hysteresis of A 6 = 25°. The raising 
of a second peak in the PDF at large values from an ini¬ 
tially monotonously decaying function can be explained 
by a growing population of bridges with receding con¬ 
tact lines. Shortly after preparation, all bridges exhibit 
a contact angle close to ( 2 a , and the average contact an¬ 
gle (9) of bridges with large volume gradually decreases 
with time, cf. Fig [5] b). Bridges with large volume and 
a capillary pressure close to zero shrink at a smaller rate 
compared to the rate at which small bridges with ad¬ 
vancing contact angles and large negative pressure grow. 
In the final state, the majority of bridges is still in the 
advancing state and make up for the large peak toward 
smaller volumes, as can be seen from Fig. [5] c). Most 
of the remaining bridges exhibit pinned or receding con¬ 
tact lines, accounting for the second, smaller peak toward 
larger volumes. 









0 0.02 y 0.04 0.06 


FIG. 5. (Color online) Distributions for data in network 
model (NW) shown in Fig. [3] a) PDF of bridge volume V 
and b) Average contact angle (9) of bridges with volume V. 
c) PDF to find bridges in bins of contact angles 9 close to the 
advancing, receding, or intermediate range at different times 
t. 


C. Large time asymptotics 



FIG. 6. (Color online) a) Squared variance of capillary pres¬ 
sure, ffp, and b) of bridge volumes, a„, as a function of time 
t for A 9 = 0° (red solid lines) and A 9 = 25° in the network 
model and in the mean field approximation (dashed lines) for 
setup B. Insets show same data at early times in linear scale. 


ment with the dimensionless cross-over time observed in 
experiments and in set up A. At early times t <C r e the 
evolution of both a p and cr^ observed in the full network 
model is well captured by the mean field model. This 
is not surprising as the capillary pressure of bridges is 
not correlated immediately after the initialization. The 
build-up of spatial correlations in the capillary pressure 
at later times t r e will invalidate the mean field model. 


Inspection of the direct comparison of the full network 
model, the mean field model and the experimental data 
of Ref. [13 in Fig. [3] demonstrates that the mean field 
model predicts the dynamics for the chosen contact angle 
hysteresis of Ad = 25° at early times rather well. In order 
to further assess the range of validity of the mean field 
model, and for a quantities comparison of the evolution 
of capillary pressures and bridge volume at late times we 
will from now on consider the larger setup B. 

Figure [S] a) and[S]b) display a comparison of the square 
of variances a p and ay corresponding to the PDF of cap¬ 
illary pressure P and bridge volumes V, respectively, in 
the mean field model (dashed lines) and the full network 
model (solid lines). In contrast to the data shown in 
Fig. [3l we use the non-dimensionalized time in Fig. [6] As 
for the simulations of set up A we chose a contact angle 
hysteresis of Ad = 0° (red solid lines) and 25° (green 
solid lines). 

In case of a vanishing contact angle hysteresis, we find a 
power-law decay of both a p and in t with an exponent 
—3/2, as expected for continuum models for linear diffu¬ 
sion in three spatial dimensions, cf. the red solid lines in 
Fig. |H]a) and[H]b) @|. The cross-over from the initial ex¬ 
ponential decay to the power-law decay can be observed 
at an estimated time t = T e ~ 10~ 3 . This is good agree- 


The decay of the squared variance a p for a contact 
angle hysteresis Ad = 25° is slightly faster than for the 
case of vanishing contact angle hysteresis (cf. the green 
and red solid line in Fig. [6] a). A pronounced qualitative 
difference, however is observed for the squared variance 
of volumes a %, displayed as the green and red solid lines 
in Fig. [G] b). 

Inspection of Fig. [G] b) shows that the volume fluctu¬ 
ations in the initial bridge volume do not fully decay in 
the limit t oo if the contact angle hysteresis is set to 
the finite value Ad = 25°. Instead the variance of bridge 
volumes shown by the green solid curve in Fig. [6]b) sat¬ 
urates to a value cr p (t = oo) > 0. Apparently, the time 
at the cross-over and the level of the plateau is reason¬ 
ably well captured by the mean held model (dashed black 
line). Here, we would expect to observe a break-down 
of the mean held model at late times because the cap¬ 
illary pressure fluctuations become correlated irrespec¬ 
tive of the magnitude of contact angle hysteresis. The 
history dependence of the contact angle, however, pre¬ 
vents a build-up of large correlations in the fluctuations 
of neighboring bridge volumes. The time of the cross¬ 
over of (Ty to the plateau is comparable to r e obtained 
from simulations for Ad = 0°. 
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D. Cross over time 


The cross-over time r e for the case without contact 
angle hysteresis can be related to the time when the cap¬ 
illary pressure of neighboring bridges starts to be cor¬ 
related. To obtain an estimate for r e we employ the 
continuum picture of a mode expansion of volume fluc¬ 
tuations into a spectrum of modes with short and long 
wavelengths. Similar arguments have been applied to 
determine the cross-over time of diffusive solvent trans¬ 
port in a network of densely packed emulsion droplets 
[271 ]. Provided the equilibration dynamics is described 
by a linear equation, we can think of the spatial fluctu¬ 
ations of bridge volumes as a superposition of modes in 
the local liquid saturation c. According to the Nyquist- 
Shannon sampling theorem [28|, the shortest meaningful 
wavelength of such a mode c(r, t) is given by twice the 
smallest distance of two neighboring bridges. Toward 
large wavelengths, the spectrum is limited by the system 
size L. 

Solutions of the linear diffusion equation D e Ac = dtc , 
with the three dimensional Laplacian A and the effective 
diffusion constant D e , show that amplitudes c q of a mode 
c(r,t) = c q exp(iq ■ r — i/r q ) decay exponentially with 
a timescale r q = D~ x |q| -2 . The characteristic decay 
time of the mode with the largest wave number g max 
provides us with an estimate of r e . Employing 2Rq as 
the typical distance between two neighboring bridges we 
have g max ss 7t/4 R 0 . 

An estimate of the effective diffusion constant D e can 
be obtained from the link between the continuum Diffu¬ 
sion equation and the discrete network model eqn. 03- 
Considering the right hand side of eqn. m as an ap¬ 
proximation of the Laplace operator we need to correct 
the left hand side of m for a numerical prefactor given 
by the coordination of neighboring bridges for a proper 
mapping onto the continuum diffusion equation. For the 
regular bridge network of setup B we find a number of 
N c = 10 neighbors. When expanding the capillary pres¬ 
sure Pi in the volume fluctuation around the final volume 
V°° up to linear order we obtain a linear diffusion equa¬ 
tion for the fluctuations of the saturation c(r, f). 

After collecting all factors, we find that the effective 
diffusion constant D e of volume fluctuation in the bridge 
network is given by an expression 


D e 


Nc dP 

2d dV 


v=v a 


( 20 ) 


where d is the spatial dimension of the bridge network. 
Setting global parameters separation S = 0, contact an¬ 
gle 6 = 7°, asymptotic volume V°° ss 0.02 and the con¬ 
ductance coefficient (7=1, we obtain an estimated ef¬ 
fective diffusion constant D e ss 500 and a corresponding 
timescale r e « 10~ 3 which compares well with the ob¬ 
served cross-over time in simulation of set-up B, cf. the 
beginning of the section IIII Al We find a faster diffu¬ 
sive spreading in the network with a diffusion constant 
D e ss 10 3 if we consider capillary bridges with a pinned 



FIG. 7. (Color online) Temporal evolution of the squared 
variance cr 2 of bridge volumes in the full network model 
for different magnitude of contact angle hysteresis Ad = 
0°, 4°, 8°, 15°, 25°. Inset: Double logarithmic plot of cr 2 (t = 
oo) against Ad in comparison to a power law ~ Ad 3 (dashed 
line). 


contact line and opening angle 0 = 21.5° (corresponding 
to a volume of U 00 = 0.02). We obtain this estimate by 
replacing the capillary pressure P in eqn. for bridges 
with freely sliding contact lines at a fixed contact an¬ 
gle with the pressure P for bridges with pinned contact 
lines. An enlarged value of the derivatives of the capil¬ 
lary pressure with respect the volume is also apparent in 
the sketch in Fig. [2] The accelerated diffusive transport 
of liquid between capillary bridges with pinned contact 
lines becomes even more pronounced for larger contact 
angles. 

So far, we have explored the differences in the late time 
regimes of equilibration of bridge volume for the two ex¬ 
treme cases of a large contact angle hysteresis Ad = 25° 
and vanishing hysteresis. A cross-over to a power-law 
scaling of the variance of bridge volumes observed for the 
ideal case Ad = 0° can be understood in the framework 
of a continuum model for diffusion. To quantify the in¬ 
termediate cases we considered the late time regimes for 
a series of contact angle hysteresis. Figure [7] shows the 
decay of the square of the volume variance, <j 2 , for addi¬ 
tional values Ad = 4°, 8°, 15° in set up B. The plot in 
the inset of Fig. [7] displays the value <r 2 (t = oo) at sat¬ 
uration in a double logarithmic plot. The available data 
suggest a power law of the cr 2 (f = oo) oc Ad 3 that needs 
to be confirmed in further investigations including larger 
regions of the parameter space. 


IV. OUTLOOK AND CONCLUSIONS 


The present numerical study demonstrates that con¬ 
tact angle hysteresis strongly affects the evolution and 
distribution of capillary bridge volumes after mixing a 
wetting fluid into a bed of spherical beads. Our model 
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study shows that disorder of the local geometry in the 
granular bed such as the distribution of gap separations 
and the polydispersity of bead radii hardly contribute to 
the final width of bridge volume distribution. Unless the 
difference between the advancing and receding contact 
angle on the beads is unrealistically small, contact angle 
hysteresis will always dominate the final distribution of 
bridge volumes. 

Our numerical results further demonstrate that a mean 
field model for the capillary pressure is sufficient to de¬ 
scribe the equilibration process at early times. A reli¬ 
able prediction of the final distribution of bridge volumes, 
however, must account for the network of communicat¬ 
ing bridges. The time of the cross-over between the early 
time regime with negligible spatial correlations of the vol¬ 
ume fluctuations to a late time regime where a continuum 
model describes the volume distribution can be derived 
for the ideal case of a vanishing contact angle hysteresis. 
This cross-over time provides an estimate for the satura¬ 
tion time of spatial volume fluctuations in case of high 
contact angle hysteresis. 

The network model can be further employed to address 
the role of contact angle hysteresis in fluid transport in 
slowly sheared granular beds, or the dynamics of liquid 
equilibration in particulate materials driven by gradients 
in the local saturation. It can be easily extended to ac¬ 
count for flows of the wetting fluid induced by gravity or 
evaporation. 


V. APPENDIX 

Starting point for the derivation of eqns. (Q} and (0 for 
a one component system is the condition of coexistence 

Hv(p v ,T) = fii{pi,T) (21) 


an ideal gas, we can employ the relation p v /p* = n v /n* 
and rewrite eqn. (12411 in the form 


P n. 


*2 




p* nj 


(25) 


in terms of the capillary pressure P = pi — p v . Equation 
m is equivalent to Kelvin’s equation. 

In most cases, the liquid coexists with a mixture of the 
liquid vapor and an inert gas that is insoluble in the liquid 
phase. In this case we need to replace the equilibrium 
condition eqn. (EU) by 


Pv(Pv,T) = pi(p 0 + P,T) (26) 


where po is the ambient atmospheric pressure, p v the 
partial pressure of the liquid vapor. An expansion of 
both sides of eqn. m around the partial vapor pressure 
p v of a flat interface defined by 

Vv(pv,T) = m(po,T) (27) 

provides us with the Kelvin equation 

Pn 2 

n v -n v = ^ , (28) 

PvTll 

which is of the same form as eqn. (HT51) . All bulk quan¬ 
tity occuring in eqn. (BHD , however, refer to a constant 
ambient pressure po of the gas phase. 

According to Fick’s first law of diffusion, the current 
j of liquid molecules in the vapor phase can be related 
to the concentration gradient V n v and the diffusion con¬ 
stant D by 


j = -D V n v . (29) 


for the chemical potential /q of bulk phases i G {u,/} 
at temperature T. Due to the curved interface between 
the liquid and vapor phase, we have pi ^ p v ■ Expanding 
the chemical potential /q(p, T) of the liquid and vapor 
phase around the pressure p * at coexistence with a flat 
interface leads to 

A m(pi,T) » p*(T) + ^ Api . (22) 

d P P = P * 

with A pi = Pi — p*. 

Maxwell’s relation for the grand potential 
Gi(pi,N,T) = NiPi{pi,T) of the bulk phases i G {v,l} 
provides us with the expression 


dpi _ l_ 
dp Hi 


(23) 


for the molar density rq of the bulk phases, and we can 
rewrite eqn. d in the form 


nf — n* . nf . 

Pi-Pv = --— A p v « A p v , (24) 

K K 

where the last approximation is justified as long as 
nf /nf 1 holds. Since we assume the dilute vapor to be 


Since the total particle current Njj between bridge i and 
j scales with the interfacial area of the bridges ~ Rq 
and with the gradient |V| ~ -RcT 1 ; we finally arrive a at 
volume flux of the liquid phase 

Q>, = Q, Dr f*° (P~Pi) > (30) 

p* n* 

where Rq is the bead radius and Cq- a numerical prefactor 
that accounts for the particular geometry of the bridges 
and pore space. 
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